Lyapunov instability of rough hard-disk fluids 
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The dynamical instability of rough hard-disk fluids in two dimensions is characterized through 
the Lyapunov spectrum and the Kolmogorov-Sinai entropy, /iks, for a wide range of densities and 
moments of inertia I. For small I the spectrum separates into translation-dominated and rotation- 
dominated parts. With increasing I the rotation-dominated part is gradually filled in at the expense 
of translation, until such a separation becomes meaningless. At any density, the rate of phase-space 
mixing, given by hns, becomes less and less effective the more the rotation affects the dynamics. 
Ifowever, the degree of dynamical chaos, measured by the maximum Lyapunov exponent, is only 
enhanced by the rotational degrees of freedom for high-density gases, but is diminished for lower 
densities. Surprisingly, no traces of Lyapunov modes were found in the spectrum for larger moments 
of inertia. The spatial locahzation of the perturbation vector associated with the maximum exponent 
however persists for any /. 



I. INTRODUCTION 



The phase space trajectory of a many-body system 
with convex particles is Lyapunov unstable [l], [2| , which 
means that it is extremely sensitive to small (infinitesi- 
mal) perturbations of the initial conditions. If all possi- 
ble perturbation vectors are represented as an infinitesi- 
mal hypersphere centered on and co-moving with a phase 
point, at later times this object is deformed under the ac- 
tion of the linearized flow. The main axes exponentially 
grow or shrink with time. The time-averaged rates are 
referred to as Lyapunov exponents, and the sorted set 
of exponents, {A;,^ = as Lyapunov spectrum. 

Here, D is the phase-space dimension. The relation of the 
Lyapunov spectrum to more familiar fluid properties has 
motivated much work over the last two decades. If, for 
example, the density of a soft-potential Weeks-Chandler- 
Anderson fluid [sl is isothermally increased from the fluid 
to the solid phase, the Kolmogorov-Sinai (dynamical) en- 
tropy, which is equal to the sum of all positive Lyapunov 
exponents exhibits a maximum at a density which is 
about 20 % smaller than the liquid-line density at the 
fluid-to-solid phase transition. This result, which has 
been verified for two- [1, Hj and three-dimensional 0, [1| 
systems, is somewhat surprising. It implies that in soft- 
potential systems most efficient phase space mixing and 
information generation about the initial state does not 
occur at the transition itself but at a slightly smaller den- 
sity. This sheds new light on our understanding of the 
onset of phase transitions. Another example concerns 
stationary non-equilibrium systems with time-reversible 
equations of motion. In this case, the Lyapunov spectrum 
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may be directly linked with the transport coefficients and 
the rate of entropy production P, . 

Systematic studies of simple atomic fluids, such as 
hard-sphere and soft-sphere models, have revealed some 
interesting features. In particular, trajectory pertur- 
bations related to the largest (in absolute value) Lya- 
punov exponents are strongly localized in physical space 
[111 . [T^ . ITsj . This means that at any instant of time only 
a very small fraction of all particles contributes to the 
fastest dynamical events responsible for the largest rates 
of perturbation growth (or decay) . This localization also 
persists in the thermodynamic limit. A similar localiza- 
tion property for the maximum-exponent perturbations 
was found for one-dimensional distributed systems with 
space-time chaos On the other hand, perturbation 
vectors associated with the smallest (in absolute value) 
exponents may be delocalized and exhibit wave-like pat- 
terns in space. These so-called Lyapunov modes were 
flrst observed for systems of hard dumbbells [H, and 
hard disks [13, [H, [II] and are a consequence of the ba- 
sic symmetries, namely invariance with respect to time 
and space translations (depending on the boundary con- 
ditions). Due to exponent degeneracies, they are recog- 
nized by a step-like appearance of the Lyapunov spec- 
trum for small (in absolute value) exponents. For soft- 
potential systems, however, sophisticated Fourier trans- 
form techniques are required to demonstrate the presence 
of Lyapunov modes [1, [l^] . 

The study of Lyapunov exponents has also been ex- 
tended to simple molecular systems such as soft [2lj, [13] 
and hard dumbbells [isl [l6j in two dimensions. In this 
case the dynamics is affected by qualitatively different 
degrees of freedom, translation and rotation, which both 
have a pronounced effect on the spectrum. 

In this paper we turn to a related model, which also 
incorporates the idea of the coupling of translational 
degrees of freedom to some internal energy, which af- 
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fects the dynamics, namely rough hard disks in two di- 
mensions. The three-dimensional version of this model, 
rough hard spheres, was first introduced by Bryan 23 1 
in 1894, and was consecutively treated by Pidduck 24 1 
and, in particular, by Chapman and Cowling |25j . who 
worked out explicit formulas for the transport coefficients 
in terms of kinetic theory. It is an extension of the famil- 
iar smooth hard-sphere model, with angular momentum 
added to each sphere due to roughness. Roughness is 
introduced by the requirement that a collision between 
two particles reverses the relative surface velocity at the 
point of contact of the collision partners, exchanging both 
linear and angular momentum in the process. This def- 
inition corresponds to the maximum possible roughness 
between two particles. This model was extensively stud- 
ied by Berne and co-workers with respect to various mod- 
els of molecular rotational relaxation [2^ . They also in- 
vestigated the dynamical properties of models with par- 
tial roughness in between that of the smooth and the 
maximum rough sphere models p7l . f28j , in particular the 
hydrodynarnic long-time behavior of various correlation 
functions ,29]. 

Here we are concerned with the two-dimensional ver- 
sion of this model with maximum roughness, N identical 
rough hard disk in a box with periodic boundaries. The 
paper is organized as follows: In Section |TT] we describe 
the model and derive the collision map for colliding par- 
ticles and the respective linearized map for the dynamics 
in tangent space. All simulation results are presented in 
Section IIIII We conclude with a discussion of the results 
in Sec. Ml 



II. 



ferential equations. 



ROUGH HARD SPHERE AND HARD DISK 
MODELS 

Phase space dynamics 



We consider the three-dimensional {d — 3) rough hard- 
sphere model first. It consists of N identical hard spheres 
of diameter a, mass m and moment of inertia /, which 
are located at positions qi, and which move with (lin- 
ear) velocities Vi and rotate with angular velocities cui, 
i e {1, • • • , N}. Due to the isotropy of the spheres, the 
particle orientation is not required in the following, and 
the state vector is given by 



f=[m,{v.},{Lj,}). 



(1) 



The phase space dynamics is characterized by force-free 
" streaming" , interrupted by instantaneous pairwise col- 
lisions, at which momentum and angular momentum are 
transferred between the colliding particles, such that the 
relative surface velocity g at the point of contact of the 
two particles is reversed. Their positions are not affected 
by the collisions. 

The equations of motion for the streaming between 
binary collisions is written as a system of first-order dif- 



r = F(f): ({gi = i/,},{^, = 0},{ti;, = 0}), (2) 

where i = 1, . . . , N . It has an obvious solution. 

At a collision between two particles i and j, the colli- 
sion map 



f ' = M(f) 



(3) 



is introduced, which relates the pre-coUision state F to 
the post-collision state F'. With the following defini- 
tions, 

Q = q]-qu n=-q\ V = Vj - vf, rt = ujj+u;i, (4) 
a 

the relative surface velocity at the point of contact of i 
and j becomes 



^ 2 



(5) 



Here, n denotes a unit vector along the line of centers 
from i to j at the instant of collision, and v is the re- 
spective relative velocity. Maximum roughness requires 
that 
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-9, 



(6) 



where the prime, here and below, refers to the state im- 
mediately after the collision. Together with the conser- 
vation laws for energy, linear momentum and angular 
momentum, that part of the collision map relevant for 
the collision between particles i and j becomes [Isj 



Vi' = Vi+-fg + (3 n{n ■ t/), 

Vj' = Vj -jg- Pn(n-v), (7) 

LUi ' = LUi + (2/?/(t) rt X g, 

ujj ' = LUj + (2(3/ a) n X g. 

Here, /3 and 7 are constants, which depend on the mo- 
ment of inertia / around a diameter of the sphere: 



4/ 



1' 



/3 = 



1 



1 



(8) 



The dimensionless parameter 7 controls the coupling be- 
tween translational and rotational degrees of freedom. In 
the limit / ^ 0, the translational and rotational degrees 
of freedom decouple and the rough hard sphere model 
reduces to the conventional smooth hard-sphere model 
without roughness [s^l , if all angular velocity vectors are 
discarded. 
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Tangent space dynamics 

A perturbed trajectory is separated from the reference 
trajectory by an offset vector 



If the moment of inertia vanishes (7 — > and (3—^1), the 
hnearized coUision map of the smooth hard sphere fluid 
is recovered [3l| , if the angular velocity perturbations are 
discarded. 



ST^i{Sq,},{Sv,},{6Lu,}), 



(9) 



which evolves during the streaming phase according to 
the linearized equations of motion 

^ dP ^ / ^ - ^ ^ 
Sr = ^-Sr: {% = 6v,}, {Sv, = 0}, {Suj^ = 0} 

dr ^ 

(10) 

i = 1, ■ ■ ■ , N. It is trivially solved. 

The linearization of the collision map ^ is obtained 
fromRef. ,301: 



dr 



^.F(f)-i?(M(f)) 



St,. (11) 



Here, 6tc is the (infinitesimal) time shift between the 
collision of the reference trajectory and of the perturbed 
trajectory, which may be positive or negative. Similarly, 
we denote by Sq, the shift in configuration space of the 
collision points of these two trajectories. If, as before, 
i and j are the colliding particles, these quantities are 
computed from [31 1 



Stc 



Sq- n 



Sqc = Sq + V Stc 



(12) 



where we use a notation for the perturbed quantities 
which is analogous to the previous notation in phase 
space (see Eq. Q): 

5q = 5qj—5qi] 5v~5vj — 5vi] 5Vt ~ Sujj + SCoi. (13) 

Since 5q, = aSn, we also find from Eq. 



dq = Sv -\ — 
2 



5qc X n + qx sn 



(14) 



With this notation, we obtain for that part of the lin- 
earized map (1111) belonging to the collision of i and j: 



Sq/ = Sq,- ig+^q{q-v) 

L (7^ 



5qj ' = 5qj + 



/9 
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q{q- v) 



Stc , 



Stc 



Svi ' — Svi + jSg 

9 



Sqdq- v) + q{v ■ Sqc) + q{q ■ 6v) 



Svj ' — Svj — jSg 



f3 



HM- v) + q{v ■ 5qc) + q{q- Sv) 



5uj 



21 
(72 

2/3 



' ' (72 



6qc X g + q X Sg 
ScfcX g + qx Sg 



(15) 



Computer simulations of hard disks in two dimensions 

For the remainder of this work we restrict ourselves to 
the case of planar rough disks on the xy plane, d ^ 2. 
All the equations above remain valid in this case, if all 
position and velocity vectors are placed in the xy plane, 
and all angular velocity vectors LOi are perpendicular to 
this plane with a single non- vanishing z component LUi. 
Discarding all superfluous components, we are left with 
D = 5N components for the state vector T and for any 
perturbation vector ST. 

Depending on the mass distribution of the disks and, 
hence, their moment of inertia / with respect to a perpen- 
dicular axis through the center, the coupling parameter 
K may take values between zero and one: k = applies, if 
all the mass is located in the center, k = 1/2 corresponds 
to a uniform mass distribution, and k = 1 is obtained, if 
all the mass is concentrated on the perimeter of the disk. 

For our numerical work we use reduced units, for which 
the disk diameter a, its mass to, and the Boltzmann con- 
stant k are set to unity. The time-averaged translational 
kinetic energy per particle, {K) /N — m{vi ■ Vi)/{2N) 
is taken as the unit of energy. Thus the temperature is 
one, T = 1, as in our previous work on smooth hard 
disks, which facilitates comparison [sO]. We have ascer- 
tained that the translational and rotational temperatures 
agree. The unit of time is {ma'^N/{K)y/^. Lyapunov 
exponents and the Kolmorov-Sinai entropy are measured 
in units of yJkT /ma'^. The number density is defined as 
p — Na^ /V , where N is the total number of disks, and 
V denotes the area of the simulation box with extensions 
Lx, Ly. For fluid phases we take a square box, = Ly, 
where the particles are initially put on a square lattice 
with random velocities and random angular velocities. 
For solid phases the aspect ratio Ly/L^ — V5/2 is used, 
which is compatible with the triangular close-packed lat- 
tice and does not differ much from a square. By initially 
putting the particles on a triangular lattice with random 
velocities and angular velocities, high density systems up 
to the close-packed density po = 1.1547 may be stud- 
ied in this way. Periodic boundary conditions are used 
throughout. The only free parameters are the number 
density p and the moment of inertia / (respective the 
coupling parameter k — 4/). 

To evolve the system, an event-driven algorithm simi- 
lar to that of Rapaport [32] was used. Proper care was 
taken to avoid missing glancing collisions. The conserva- 
tion of energy and linear momentum was carefully ver- 
ified. The total angular momentum is not conserved 
due to the periodic boundaries. For the computation 
of the 5A^ Lyapunov exponents, the classical algorithm 
of Benettin et al. [1^ and Shimada et al. [34|, properly 
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modified for the instantaneous elastic rough coUisions en- 
countered here, was used. Simultaneously with the refer- 
ence trajectory r{t), 5N rephca of perturbation vectors, 

(jf*" \t), I = 1, • • • , 5A^, with orthonormal initial condi- 
tions were evolved and periodically re-orthonormalized 
with a Gram-Schmidt procedure [30., ,31j. 



III. RESULTS 

Lyapunov spectra 

In this section we discuss the Lyapunov spectra of 
rough disk fluids and solids. A full spectrum consists 
of the ordered set of Lyapunov exponents, Ai > A2 > 
■ ■ • > A^i, where D is the phase space dimension. For 
a symplectic system as in our case, the Lyapunov expo- 
nents always appear in pairs with a vanishing pair sum. 
Therefore, the Lyapunov spectrum consists of a positive 
and a negative branch and is symmetric, A; = —Xd+i-i, 
if plotted as a function of the index /. In the absence 
of magnetic fields, this so-called "conjugate pairing sym- 
metry" may simply be viewed as a consequence of the 
time-reversibility of the dynamical equations. Therefore, 
we may restrict ourselves to the positive branch of the 
spectrum, 1 < I < D/2, for which A; > 0, which consid- 
erably reduces the computational effort. Examples for 
full spectra will be given below (see the upper panel of 
Fig. [51). 

Another important property concerns the vanishing 
exponents, which are a consequence of the fundamen- 
tal continuous symmetries which leave the Lagrangian 
and, hence, the motion equations invariant. According 
to Nother's theorem, each such symmetry corresponds 
to a constant of the motion [sst and, in addition, gen- 
erates two symplectic-conjugate vector fields in phase 
space, along which perturbations do not stretch or shrink 
and, therefore, give rise to two vanishing exponents 0. 
For example, invariance with respect to time translation 
gives rise to energy conservation and to two vanishing 
exponents. Similarly, invariance with respect to uniform 
translation in space gives rise to momentum conservation 
and four more vanishing exponents. However, isotropy of 
space and, hence, angular momentum conservation ap- 
plies locally for each collision, but not globally due to 
the periodic boundary conditions. Altogether only 6 Lya- 
punov exponents vanish in our case, as is also confirmed 
by the simulations. 

Although a spectrum only consists of discrete points in- 
dexed by the integer Z - or by the reduced index I = 21 /D 
- most of the time in the figures below it is represented 
by a smooth line drawn through these points to enhance 
the clarity. 

To assess the influence of translation-rotation coupling 
on the Lyapunov spectra, we show in Fig. [T] results for a 
rather dilute gas, p = 0.1, of iV = 400 rough disks at a 
temperature T = 1. The various curves belong to differ- 



K = 0.004 

K = 0.100 
K = 0.200 

K = 0.400 

K= 1.000 




0.16 
0.14 
0.12 
0.1 
0.08 
0.06 
0.04 
0.02 


-0.02 



K = 0.004 
K = 0.100 
K = 0.200 
K = 0.400 
K= 1.000 



750 



760 



770 



780 
/ 



790 



800 



FIG. 1: (Color online) Lyapunov spectra for a system of 
A*' = 400 rough hard disks at a low density p = 0.1. The 
various curves are for different moments of inertia /, where 
the keys refer to the coupling parameter k = AI. The uniform 
mass distribution corresponds to /t = 0.5. Top panel: Posi- 
tive branches of the specta. The reduced index I is used on 
the abscissa. Bottom: Magnification of the transition region 
between translation and rotation dominated exponents. The 
un-normalized index / is used on the abscissa. 



ent moments of inertia I and are specified by their cou- 
pling parameters k = 47. The system is fairly large, and 
the results are close to the thermodynamic limit [30| . The 
phase space has D — 2000 dimensions. In the top panel 
of Fig. [1] the positive branches of the full spectra are 
shown, where the normalized index I = I /{D/2) is used 
{\ <l < D/2 — 1000) on the abscissa. Most noticeable is 
the transition region near Iq = {2N - 'i)/{D/2) = 0.797, 
which separates the specta into a translation-dominated 
regime for / < Iq and a rotation-dominated regime for 
^0 < ^ < (5iV — 6)/D). A magnification of the transi- 
tion region is shown in the lower panel of Fig. [U where 
the un-normalized index I is used on the horizontal axis. 
Analogous spectra for a rather dense gas, p = 0.7, are 
shown in Fig. [2l 

For very small k translation and rotation are effectively 
decoupled and the dynamics is almost identical to that 
of a smooth hard-disk gas. For 400 particles the positive 
exponents Ai, • • • , A797 agree with the positive exponents 
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FIG. 2: (Color online) Lyapunov spectra for a system of 
— 400 rough hard disks at a moderately high density 
p = 0.7 (top). The various curves are for different moments of 
inertia, where the keys specify the coupling parameter k = 41. 
The uniform mass distribution corresponds to k = 0.5. Top 
panel: Positive branches of the specta. The reduced index I 
is used on the abscissa. Bottom panel: Magnification of the 
transition region between translation and rotation dominated 
exponents. The un-normalized index I is used on the abscissa. 
Reduced units and periodic boundaries are used as explained 
in the main text. 



of the smooth hard disks and, hence, are translation dom- 
inated, whereas the (very small but positive) exponents 
A798,---A997 are due to the angular velocity perturba- 
tions and are only present in the rough-disk case. The 
three remaining exponents, Aggg,--- ,Aiooo, vanish due 
to the conserved quantities. Note that this is only the 
positive branch of the full spectrum. 

The maximum Lyapunov exponent Ai is generally 
taken as in indicator and a measure for dynamical chaos. 
Similarly, the Kolmogorov-Sinai (or dynamical) entropy 
hxs is a measure of phase-space mixing [36]. Due to the 
exponential instability, a number of initially close phase 
points are eventually uniformly distributed over the en- 
ergy surface. The characteristic time for this mixing pro- 
cess is the mixing time I/Hkb [s^, [11]. Since, according 
to Pesin [1] , hxs is equal to the sum of the positive ex- 
ponents, it is directly accessible through the Lyapunov 



FIG. 3: (Color online) Dependence of the maximum Lya- 
punov exponent Ai (upper panel) and of the Kolmogorov- 
Sinai entropy per particle, Kks/N, (lower panel) on the cou- 
pling parameter k for various densities as indicated by the 
keys. The rough disk fluid contains N=400 particles and has 
(translational) temperature kT — 1. Reduced units are used 
throughout. 

spectrum. In the upper panel of Fig. [3] the maximum 
exponent as a function of k is shown for various densi- 
ties. An analogous plot for the KS-entropy per particle, 
fiKs/N, is provided in the lower panel of the same fig- 
ure. If K is increased, Ai decreases - weakly - for lower 
densities, p < 0.7, and increases for large densities. The 
KS-entropy always decreases with k, even for large den- 
sities. This means that mixing becomes less effective the 
more the rotational degrees of freedom affect the transla- 
tional dynamics. Similarly, chaos is slightly reduced with 
increasing k, at least for low-density gases. 

To demonstrate the density dependence, we show in 
the upper panel of Fig. [¥] Lyapunov spectra for various 
densities of a 400-particle gas of hard disks with a very 
small moment of inertia for which k = 0.004. Analogous 
spectra with a large moment of inertia corresponding to 
K — 0.4 are provided in the lower panel of the same figure. 
All exponents increase with the density, in particular Ai, 
as is shown in more detail in Fig. [S] where the maximum 
exponent - for various k - is plotted as a function of p. In 
some sense, Ai behaves similar to the potential-generated 
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FIG. 4: (Color online) Lyapunov spectra (positive branch 
only) of 400 rough hard disks with very small moment of iner- 
tia {n — 41 — 0.004, upper panel) and with a large moment of 
inertia (/t — 0.4, lower panel) for various densities as specified 
by the key. The reduced index / is used on the abscissa. 



contribution to the pressure P [4l|. For low densities, 
both (P/pkT) — 1 and Ai are proportional to the sin- 
gle particle collision frequency U2 . Fig. [5] resembles the 
respective phase diagram for the pressure. The conspicu- 
ous gap in the spectra marks the two-phase region for the 
fluid-solid transition. As mentioned in the previous sec- 
tion, the data beyond the solid line were obtained with a 
non-square simulation box (aspect ratio ^/3/2), the data 
below the fluid line with a square simulation box. This 
choice of aspect ratios merely facilitates the setting up of 
the initial conditions for the solid-state simulations and 
does not have any signiflcance for the results. The gap 
disappears completely, if Ai is plotted as a function of the 
single-particle collision frequency 1^2 (not shown), which 
is easily obtained from the simulation. This has been 
noted already for the smooth hard-disk system (k = 0) 
(30| and is also true for all k > 0. This means that the 
statistical distributions for the parameters characterizing 
the collisions (such as the impact parameter) do not no- 
ticeably differ for the disordered fluid and the coexisting 
crystal. 

Based on kinetic theory, a density expansion for Ai of 
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FIG. 5: (Color online) Density dependence of the maximum 
exponent of 400 rough hard disks for various moments of in- 
ertia corresponding to the specified values of n. 

the smooth hard disk model (k = 0) becomes 

Ai =Ai/2[-lnp- 5 + 0(1/ In p)], (16) 

where 

1/2 = 2Tr^^^ pa {kT/m) g{a) (17) 

is the single-particle collision frequency. The pair dis- 
tribution function at contact, g{a), converges to unity 
in the low-density limit. Estimates for the constants A 
and B have been computed by van Zon and van Bei- 
jeren, A = 1.473,5 = 2.48 [3i,|43|, which represent the 
numerical data well for very small densities p < 10^"^. 
Expressions such as Eq. [16] with different constants are 
expected to hold also for the rough disks when k > 0. 

Whereas all of the results so far are for systems con- 
taining 400 particles, the dependence of the KS-entropy 
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FIG. 6: (Color online) Dependence of the Kolmogorov- Sinai 
entropy Hks on the single-particle collision frequency of 
= 64 rough disks for various k as indicated by the keys. 1/2 
is obtained by dividing the experimentally determined total 
number of collisions by N/2. 
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on the single-particle collision frequency V2 is demon- 
strated in Fig. [HI for a system with only iV = 64 disks (to 
reduce the computational cost) . This number is still large 
enough to be representative for large systems. As was 
found for the maximum Lyapunov exponent, the phase 
transition (near V2 ~ 20) does not specifically show up 
in Hks when viewed as a function of V2 instead of p. If 
p approaches the close-packed density po = 1.1547, both 
Ai and Hks diverge due to the divergence of 1/2 ■ For very 
low densities and smooth hard disks (k = 0), a kinetic- 
theory based density expansion for the KS-entropy simi- 
lar to that in Eq. [TBI becomes 

hKs/N ^ A'u2 [-\np + B' + 0{p)] . (18) 

Estimates for the constants A' . B' have been obtaine d by 
van Zon et al. [HI and most recently by de Wijn f 43| , 
A' = 0.5, B' = 1.47 ± 0.11, which describe the numerical 
data well for p < 10^"^. Again we expect a similar repre- 
sentation to hold also for k > 0. Still, it seems surprising 
that the KS-entropy is reduced so much by the introduc- 
tion of an internal degree of freedom (rotation), which 
effectively acts as energy storage in between collisions. 

An interesting step-like structure is observed for very 
small K in Fig. [2] and in the upper panel of Fig. |4l These 
steps are a remnant of the degeneracy of exponents due 
to the existence of Lyapunov modes for smooth {k = 0) 
hard disk systems [13, [3 ■ Lyapunov modes are periodic 
spatial perturbations associated with the small positive 
exponents with indices I < Iq = 2N — 3 (and with the 
conjugately paired negative exponents). The correspond- 
ing perturbation vectors may be represented as periodic 
vector fields coherently spread out over the simulation 
box and with well defined wave vectors. They may be 
understood as Goldstone modes of a system with contin- 
uous symmetries (4^ - translation invariance in space and 
time - which give rise to conservation of energy and lin- 
ear momentum [35| and in addition to the six vanishing 
exponents Note that angular momentum is globally 
not conserved according to the periodic boundaries and 
does not contribute. Fig. [2] shows that the exponent de- 
generacy and, hence the Lyapunov modes are still rather 
well developed for k = 0.004 corresponding to a moment 
of inertia / — 0.001. This is independent of the den- 
sity as is shown in the upper panel of Fig. U) However, 
if K is increased, the steps quickly disappear and Lya- 
punov modes do not seem to exist any more. This is 
most clearly demonstrated in Fig. [2l It is not clear to 
us why this happens in view of the fact that modes are 
readily found for two-dimensional hard dumbbell fluids. 
Fourier transformation techniques will be required to set- 
tle this point. 

For rough hard disks the angular velocity subspace of 
the full phase space has N dimensions and contributes 
N exponents to the full spectrum. For small k these ex- 
ponents are different from zero, but small. Half of them 
belong to the positive branch (to which we restrict our- 
selves without loss of generality) and are located in the 
index interval 2iV - 2 < Z < (57V/2) - 3, sandwiched 



between the translation-dominated regime, / < 2N — 3, 
and the three vanishing exponents still attributed to the 
positive branch, (5A^/2) -2<l < 5N/2. We refer to this 
regime as rotation dominated. If k is increased, the expo- 
nents in this regime are increased, and the exponents in 
the translation dominated regime become smaller, until 
the spectrum becomes very uniform as, for example, in 
the lower panel of Fig. 21 and the separation into transla- 
tion and rotation dominated regimes becomes meaning- 
less. Such a system we call fully coupled. Translation 
and rotation contribute indistinguishably to the mixing 
process in phase space. 



Localization of tangent-space perturbations 

The maximum (minimum) Lyapunov exponent is the 
rate constant for the fastest growth (decay) of a phase- 
space perturbation and is dominated by the fastest dy- 
namical events, binary collisions. It is not too surpris- 
ing that the associated tangent vector components are 
significantly different from zero for only a few strongly- 
interacting particles at any instant of time. Thus, the 
respective perturbations are strongly localized in phys- 
ical space. It has been shown that for both hard and 
soft disk respective sphere systems the localization per- 
sists in the thermodynamic limit, such that the fraction 
of tangent-vector components contributing to the gen- 
eration of Al follows a power law oc A^~^ rj > 0, and 
converges to zero for TV ^ 00 0, [H, [H, . The local- 
ization becomes gradually worse for larger indices Z > 1 , 
until it ceases to exist and (almost) all particles collec- 
tively contribute to the coherent Lyapunov modes men- 
tioned in the previous section. Similar observations for 
spatially extended systems have been made by various 
authors 13, 18, 46, 47, 4|. 

To demonstrate the localization property of the rough 
hard-disk system, we define the contribution of an indi- 

7? (0 

vidual disk i to the perturbation vector ST belonging 

^ (0 

to A/ as the square of the projection of oT onto the 
subspace of this disk. 

Because ST is normalized in the Gram-Schmidt step of 

the algorithm, one has J^f ^^f^ = 1 ^o^' ^ ^^'^ '^^Y 
be interpreted as a kind of (normalized) action probabil- 
ity of i for the perturbation I. It should be noted that for 
the definition of pf^ the Euclidean norm is used and that 
all localication measures depend on this choice. Qualita- 
tively, this is sufficient to demonstrate localization. 

In the top panel of Fig. [71 p^P for Z = 1 is plotted 
as a scalar field p'^^\q), where the surface is interpolated 
over regular grid points covering the whole simulation 
box. There exist one big and a few smaller competing 
active zones, which move around randomly such that the 
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l/{D/2) 

FIG. 7: (Color online) Top: Instantaneous view of the action 
probability ^^^^ for the maximally localized tangent vector, 
plotted as a scalar field over the simulation box. For the grid 
interpolation a weight function ■w{r) oc is used, where 
r is the distance of a particle to the grid point. The system 
consists of TV = 144 rough hard disks with k = 0.5 at a density 
p = 0.5. 

Middle: As in the top panel, but for /i^^""' , which belongs to 
a delocalized vector with normalized index I = 0.556. 
Bottom: Localization width W^lj (see Eq. JlSJ) for iV = 144 
disks with k indicated by the labels. The density p — 0.5. 
The reduced index / = l/{D/2) is used on the abscissa, and 
D = 5iV. 

system remains homogeneous when viewed over a long 
time. This should be contrasted with the middle panel, 
where an analogous plot (with the same scale) is shown 
for the delocalized tangent vector with an index I = 200 
{I = 0.556). For this comparison, the system contains 144 
rough disks with k = 0.5 at a density p — 0.5. 

A number of localization measures have been intro- 
duced to assess the localization of ST , not only for I — 1 



[T^ but for all I [T^ [T^. The most common is due to 
Taniguchi and Morriss '1^, who define a "localization 
width" 

W^l = exp[5«]/7V, (19) 

which is based on the Shannon entropy for the "proba- 
bility" distribution 

5(0 = ^-fj^f)ln^f)y 

Here, (•••) denotes a time average. w!ji\,j is bounded 
according to 1/iV < W^lj < 1, where the lower and 
upper bounds apply for complete localization and de- 
localization, respectively. In the bottom panel of Fig. 
[7] we plot W^Ij for the 144-disk system used before. 
The value of k is indicated by the labels. This local- 
ization spectrum changes surprisingly little when k is in- 
creased from zero to one. The only major difference is 
for the rotation-dominated regime. For the smooth disks, 
AC = 0, the points for are irrelevant in this regime, 

0.8 < I < 1.2, and are not shown. Note that only data for 
the positive branch of the Lyapunov spectrum are shown, 
l< 1. 

Alternatively, an even simpler definition may be used, 
which involves the Fermi entropy (sometimes also re- 
ferred to as the quadratic entropy), [53| 

Sf^(t^'il-^^'^. (20) 

It has the desired property: Sp'^ vanishes, if only a single 
particle is responsible for the phase-space growth (ex- 
treme localization), it is (N — l)/iV « 1, if all particles 
contribute identically (complete coherent delocalization), 
and it is in between otherwise. This measure might be 
particularly useful, whenever localization is even more 
complete than in the case presented here, but it distin- 
guishes poorly between very delocahzed states.. 

At this point, a critical remark is in order. The lo- 
calization spectrum in the bottom panel of Fig. [7] is 
shown for the positive branch of the Lyapunov spetrum 
only. It should be completely symmetrical with respect 
to the conjugate negative branch, S'(') = to 
the time reversible phase-space structure: a time rever- 
sal operation converts the stable manifold into the un- 
stable manifold and vice versa. For the smooth hard- 
disk case, K = 0, this symmetry is observed with high 
numerical precision [5l|, [H^]. However, for k > the 
spectra are slightly asymmetric (not shown) . The reason 
for this asymmetry is subtle. The perturbation vectors 
-> (0 

6r we use in this work are ortho-normal. They span 
the correct subspaces of the tangent space required for 
the computation of the Lyapunov exponents according 
to the standard algorithm [l^, [H, [s^l , but they are not 
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covariant: that means, they do not strictly follow the 
linearized dynamics in tangent space, but are regularly 
re-orthonormalized by the Gram-Schmidt procedure. As 
a consequence, they are not invariant under time rever- 
sal. The last property, however, is required for a com- 
plete symmetry of the localization spectrum, such that 
- (0 . 

the expanding vector ST in the time-forward direction 

(D+l-l) 

becomes the contracting vector ST in the time- 

backward direction. If proper covariant Lyapunov vec- 
tors [5^ are used instead of the Gram-Schmidt vectors, 
the symmetry is re-established for k > 0. Details will be 
communicated in a forthcoming publication [52] . 



spheres [305 and means that for iV — > 00, at constant den- 
sity, the Lyapunov spectra quickly converge to a limiting 
curve when plotted as a function of the reduced index 
I — l/{D/2). For the rough hard disks this is demon- 
strated in the upper panel of Fig. [8l The spectra there 
are for 16 to 256 particles at a density p — 0.7, and for 
K = 0.4. Full spectra with their positive and negative 
branches are shown, which are related by the conjugate 
pairing symmetry as was mentioned in Sec. Thus, the 
derivative of the spectrum with respect to the normalized 
index exists in that limit, 

dl JV-.00 Al 



Convergence times in tangent space 

For the Lyapunov exponents to converge, the orthonor- 
mal set of tangent vectors needs to reach its proper orien- 
tation in tangent space starting from an arbitrary initial 
orientation. The convergence time varies with the num- 
ber of particles and with the index For smooth particle 
systems it was shown in Ref. [H^ that the vector asso- 
ciated with the maximum exponent aligns with a con- 
vergence time proportional to iV", where a lies between 
0.4 (for smooth hard disks) and 0.9 (for soft repulsive in- 
teraction potentials). For higher indices the convergence 
is an even slower collective phenomenon. In this section 
the methods of Ref. Q are adapted to determine the 
system-size dependence of the convergence times of rough 
disks not only for / = 1 but for all tangent vectors. 

We consider M randomly oriented orthonormal sets 

^ (0 

of tangent vectors, Am = {6r,^},l — I,-- - ,D, where 
m = 1, • • • , M. Each set spans the full D-dimensional 
tangent space and acts as an initial condition for the com- 
putation of a full Lyapunov spectrum for the same ref- 
erence trajectory. All spectra eventually converge. Any 
two tangent vectors giving rise to the same Lyapunov 
exponent but belonging to different initial sets A need 
to become parallel or anti-parallel in the course of time, 
such that their dot-product approaches ±1. To measure 
the convergence time for a given we average over all 
such possible products. 



M-l 



M 



Xi{t) = 



M{M - 1) 



E E 

m— 1 7n'— m+1 



■ 5V^, 



(21) 



Xi{i) increases with time from xa to unity, where the 
initial value (xo ~ 0.06 for N — 36) is independent of I 
due to the random orientation of the sets A and converges 
to zero for N ~^ 00. The time for which xi{t) crosses a 
threshold Q for the first time is taken as a measure of the 
convergence time r/. In the following we take O = 0.9. 
Any other choice for this threshold only results in times 
which differ by a constant factor. 

Before continuing the discussion, we note that the 
Lyapunov spectrum exists in the thermodynamic limit. 
This has been shown for smooth hard disks and hard 



^ 




0.8 I 1.2 

l = l/(5N/2) 

FIG. 8: (Color online) Top: Full Lyapunov spectra of rough 
hard disk systems for various system sizes A'^, plotted as a 
function of the normalized exponent I = l/{5N/2). All sys- 
tems have a density p — 0.7, and k = 0.4. The spectra quickly 
converge to a smooth limit spectrum for A'^ — > 00. Bottom: 
Lyapunov vector convergence times, divided by A*', for two 
system sizes, A'' = 36 and A^ = 72 as indicated by the keys. 
The points are direct simulation results, the smooth lines were 
computed from the Lyapunov spectrum via Eq. (|23|l . where 
A = 2.85. 

It has been argued in Ref. [s^l that for k = the 
decay time for the correlation function Xi{t) concerning 
the maximum exponent is determined by 1/ (Ai — A2) and, 
hence, by the inverse "slope" of the spectrum at ^ = 1. 
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These arguments also apply to all the other exponents 
such that one expects 



A 



l+l 



X, 



(23) 



to hold, where A is a fitting parameter which, for the 
choice 9 — 0.9, becomes 2.85. Rewriting this expression 
in terms of the reduced index 1 = l/ihN/2) gives 



Tl 



^NA 



AA(/)/A/ A? 



d\{l) 



dl 



(24) 



where we have replaced the finite differences by the re- 
spective differentials. Since the limiting spectral slope, 
Eq. ([m, is independent of N, the convergence time for 
any I is expected to be proportional to the particle num- 
ber N. 

Our results for ti/N are depicted in the lower panel of 
FigurelH where experimental results for A'^ = 36 and N = 
72 rough disks are shown by the points. Both systems 
have a density p — 0.7, and k = 0.4. Clearly the points 
for different N collapse onto a single curve proving the 
proportionality of r/ to N. If r/ is computed from the 
slope of the Lyapunov spectrum according to Eqs. P5|) 
or ([24]) . the smooth lines are obtained. Their agreement 
with the simulation points supports our assertion in Eq. 
(ESI). 

The symmetry obviously exhibited by the convergence 
time, Tl = td+1-i, is surprising in view of the fact that 
the algorithm treats successive vectors successively: The 
orientation of the second vector is affected by that of the 
first, that of the third vector by that of the first and 
second, and so on. 



IV. CONCLUSIONS 

In this paper we investigate rough hard-disk systems, 
arguably the most simple models of a molecular fluid with 
translational and rotational degrees of freedom. The ro- 
tation of the particles may be viewed as a mechanism 
to store internal energy, which is returned to the trans- 
lational degrees of motion with some delay. We com- 
pute Lyapunov spectra and study the effect of rotation- 
translation coupling on the dynamical stability of such 
systems. 

If the moment of inertia / of the disks vanishes, the 
translational dynamics is completely decoupled from the 
rotational degrees of freedom and the results for the 
smooth hard-disk system are reproduced. If / respec- 
tive the more relevant coupling parameter k = 4/ is in- 
creased, the Lyapunov spectrum changes drastically with 
the rotation-dominated parts of the spectrum being grad- 
ually filled in, until a separation into translation- and 
rotation-dominated parts becomes meaningless. 

The maximum exponent, Ai, which is taken as an in- 
dicator for dynamical chaos, increases with increasing k 



for large enough densities {p > 0.7), but decreases for 
smaller densities. At the same time, the Kolmogorov- 
Sinai entropy hxs always decreases. The latter, which is 
the sum of all positive exponents, gives the rate of mix- 
ing in phase space, which becomes less and less effective 
the more important the rotation is for the dynamics. We 
encounter the unexpected situation that for large den- 
sities dynamical chaos may increases with k whereas at 
the same time phase-space mixing takes longer. This 
should be contrasted to the behavior of a system of hard 
dumbbells [lB| . For a uniform mass distribution of the 
dumbbells (corresponding to k = 0.5 for the rough disks), 
both Ai and Hks increase with the molecular anisotropy, 
and the mixing time decreases. From this point of view, 
the rough disk model seems artificial. 

Another surprise is the seeming lack of Lyapunov 
modes for the rough disks with non-vanishing k, given 
the fact that modes were readily found for hard-dumbbell 
systems jl^. As for soft interaction potential systems, 
Fourier transformation methods may still give evidence 
for modes. This point deserves further investigation. 
However, the localization in physical space of the pertur- 
bation vectors associated with the maximum exponent is 
as expected. 

The localization spectrum shown in the bottom panel 
of Fig. [7] is an application of a projection of the tangent 
vectors onto the phase space of individual disks. Due to 
the time- reversal symmetry of the evolution equations, 
such projections should show definite symmetries with 
respect to the positive and negative (not included in Fig. 
[7]) branches of the Lyapunov spectrum. However, more 
often than not, these symmetries are numerically not re- 
covered by the classical algorithm. The explanation lies 
in the fact that Gram-Schmidt-orthonormalized tangent 
vectors span the proper subspaces for the computation 
of the exponents, but are not covariant with the tangent 
flow [5^. If covariant vectors are used, these spurious 
asymmetries disappear |52|. 

Up to now little is known about the mechanism gov- 
erning the convergence of tangent vectors towards their 
proper directions. For the rough hard disk systems it is 
shown that the convergence time for all I vary linearly 
with the system size, N. Furthermore, they are related 
to the slope of the spectrum at a particular I. This view 
is suggested by the existence of the thermodynamic limit 
of the specrum. 
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